Insights into phylogenetic positions and distribution patterns: Complete mitogenomes of two sympatric Asian horned toads in Boulenophrys (Anura: Megophryidae)

Abstract Boulenophrys sangzhiensis and Boulenophrys tuberogranulata, two narrow‐distributed toad species within the Megophryidae family in southern China, are experiencing population declines due to habitat loss and degradation. Despite their critical conservation status, the two species remain largely overlooked in public and scientific spheres. This study presented the first sequencing, assembly, and annotation of the complete mitogenomes of both species using next‐generation sequencing. The mitogenome of B. sangzhiensis was 16,950 bp, while that of B. tuberogranulata was 16,841 bp, each comprising 13 protein‐coding genes (PCGs), 22 transfer RNA genes (tRNAs), two ribosomal RNA genes (rRNAs), and a noncoding control region (D‐loop). The gene content, nucleotide composition, and evolutionary rates of each mitogenome were analyzed. Both mitogenomes exhibited negative AT skew and GC skew with high A + T content. ATP8 exhibited the highest evolutionary rate, while COI had the lowest. A phylogenetic analysis based on 28 mitogenomes revealed two major clades of Megophryidae, supporting the classification of two subfamilies, Megophryinae and Leptobrachiinae. Within the subfamily Megophryinae, the genus Boulenophrys was divided into two species groups. Intriguingly, despite coexisting in Zhangjiajie City, B. sangzhiensis and B. tuberogranulata exhibited distinct origins from the two different species groups, underscoring the unique role of the coexisting area Zhangjiajie in driving their speciation and preserving their current populations. A parallel pattern was also identified in the Leptobrachiinae genus Leptobrachium within the same region. This study provided valuable data references and enhanced our understanding of the molecular characteristics of these threatened amphibian species.


| INTRODUC TI ON
The toad family Megophryidae (Bonaparte, 1850) is an Asian endemic group that primarily distributed from eastern Pakistan to western China, and southward to the Philippines and the Greater Sunda Islands (Frost, 2023).Although the taxonomic validity of Megophryidae has been challenged for a long time in the last century (Zhou et al., 2023), it is now generally recognized as a monophyletic group and includes two subfamilies, Megophryinae and Leptobrachiinae (AmphibiaWeb, 2023;Frost, 2023).According to the latest statistics, the two subfamilies comprise 134 and 187 species, respectively (Frost, 2023), making Megophryidae one of the most diverse groups of amphibians in Asia as well as in the world.
The subfamily Megophryinae (Bonaparte, 1850), also known as Asian horned (or spadefoot) toads, basically reflects the distribution range and pattern of the whole Megophryidae family (Frost, 2023).
Due to significant morphological diversity and similarity, especially the lack of easily recognizable generic diagnostic characters, debates regarding generic classifications within Megophryinae have persisted for decades.Several different generic hypotheses have been proposed from various studies, including five-genus classifications (Chen et al., 2017;Delorme et al., 2006), seven-genus classifications (Dubois et al., 2021;Fei & Ye, 2016;Lyu et al., 2021;Qi et al., 2021), or a single genus proposal encompassing all Megophryinae species together as a large group (Mahony et al., 2017).These specific generic classifications have attracted dedicated followings in later studies, being embraced or revised by databases over time.Among them, the proposal of a single large genus by Mahony et al. (2017) has been widely accepted in various studies (e.g., Gao et al., 2022;Liu et al., 2018;Tapley et al., 2021;Wang et al., 2019Wang et al., , 2020)).However, a recent study has introduced a 10-genus classification from an integrative taxonomic perspective (Lyu et al., 2023).This new classification is revealed by 10 well-supported molecular phylogenetic clades and presented with combined morphological diagnoses for each genus, representing the most reliable and the latest taxonomic recognition of Megophryinae.These genera include Atympanophrys, Boulenophrys, Brachytarsophrys, Grillitschia, Jingophrys, Megophrys, Ophryophryne, Pelobatrachus, Sarawakiphrys, and Xenophrys.This 10-genus classification of Megophryinae has been adopted shortly in mainstream databases such as Amphibian Species of the World (https:// amphi bians ofthe world.amnh.org/ ) and AmphibiaChina (https:// www.amphi biach ina.org/ ), which is also the updated taxonomic system we used in this study.
The genus Boulenophrys (Fei & Ye, 2016), which includes half of the recognized species within Megophryinae (67 out of 134 total species), is the largest genus in Megophryinae (Frost, 2023).China hosts the most diverse population of Megophryinae, housing over two-thirds of the total cataloged species (Lyu et al., 2023).Unsurprisingly, species in Boulenophrys are predominantly found in southern China (61 species, see Figure 1), with a few extension westwards into Myanmar and southwards into northernmost Indochina, including Vietnam, Laos, and Thailand (Lyu et al., 2023).
Similarly, Boulenophrys tuberogranulata (Mo et al., 2010) is restricted to Tianzishan Nature Reserve and Mt.Tianping in northwestern Hunan Province, found at elevations ranging from 1000 to 1380 m a.s.l.(Lyu et al., 2023;Mo et al., 2010).Notably, B. sangzhiensis and B. tuberogranulata share overlapping distributions in Mt.Tianping.Despite this, they exhibit distinct morphological characteristics, suggesting a potentially early derivation from respective evolutionary lineages.Both species, however, face population declines due to threats such as habitat loss or degradation (Gao et al., 2022).The IUCN SSC Amphibian Specialist Group assessed their status in 2020, classifying B. sangzhiensis as Critically Endangered (CR) and B. tuberogranulata as Endangered (EN) on the IUCN Red List of Threatened Species (IUCN SSC Amphibian Specialist Group, 2020a, 2020b).
As endemic yet threatened amphibians, both B. sangzhiensis and B. tuberogranulata have received very limited public attentions and scientific concerns.Knowledges about these species in science have primarily been limited to initial morphological descriptions and distributions, with little engagement in phylogenetic studies, let alone investigations of their population structure and genetic diversity.Moreover, given the existence of over 60 species in Boulenophrys, constructing a robust phylogenetic tree using mitogenomes, for example, could significantly enhance our understanding of the phylogeny, trait evolution, and current distribution patterns within this group.However, the available mitogenome data remain quite limited.In this study, we utilized next-generation sequencing (NGS) technology to sequence and assemble the complete mitogenomes of B. sangzhiensis and B. tuberogranulata for the first time.
Our primary aim was to conduct a detailed analysis and description of these two complete mitogenomes to enhance the mitogenomic data of Boulenophrys.Then, we reconstructed a phylogenetic tree of Megophryidae based on all available mitogenomes in GenBank to scrutinize the current understanding of the phylogeny of this group.This study is expected to provide valuable data references for population genetics and conservation biology studies on these two threatened species.It will also contribute to the phylogenetic knowledge of Boulenophrys and other Asian horned toads in the future.

| Sample collection and sequencing
Permissions for the field survey in this study were obtained for scientific purposes from the local administrations, and the sample col-

| Sequence assembly, annotation, and analysis
The complete mitogenomes of B. sangzhiensis and B. tuberogranulata were assembled using NOVOPlasty 4.3 (Dierckxsens et al., 2016), with the mitogenome of Boulenophrys jingganggensis serving as a reference.The assembled sequences were then queried against the standard database of NCBI using the BLAST online program to identify highly similar sequences.The positions and directions of protein-coding genes (PCGs), ribosomal RNA genes (rRNAs), transfer RNA genes (tRNAs), and the control region (D-loop) were annotated using the MITOS Web server (Bernt et al., 2013).This annotation was cross-verified by comparing it with the information available for congeners.Further annotation and gene map visualization were carried out using the web application GeSeq (Tillich et al., 2017).

| Phylogenetic analysis
All available mitogenomes of Megophryidae species from GenBank were downloaded, and one representative sequence was selected for each species.A comprehensive dataset comprising a total of 28 Megophryidae species, including the newly sequenced B. sangzhiensis and B. tuberogranulata, was compiled for phylogenetic reconstruction.As an outgroup, Microhyla fissipes from the family Microhylidae was included.PhyloSuite (Zhang et al., 2020) was employed to extract each of the 13 PCGs from the 29 mitogenomes.Subsequently, individual PCG datasets were aligned using the MUSCLE module in MEGA X and manually checked for accuracy.Finally, these aligned and trimmed datasets were concatenated to create a combined PCGs dataset using PhyloSuite.Phylogenetic trees were reconstructed using both the maximum likelihood (ML) method with RAxML 8.0.2 (Stamatakis, 2006) and Bayesian inference (BI) method using MrBayes 3.2.7 (Ronquist & Huelsenbeck, 2003).The best-fit partitioning scheme and nucleotide substitution models for ML and BI reconstructions were determined using PartitionFinder 2 (Lanfear et al., 2017).Statistical confidence was assessed by conducting a bootstrap test with 1000 replicates for ML trees.For BI trees, posterior probabilities were calculated under a simultaneous run of 1.0 × 10 7 million generations, with sampling every 1000 generations and discarding the initial 25% of generations as burn-in.

| Mitogenome assembly, annotation, and nucleotide composition
The complete mitogenome of B. sangzhiensis was 16,950 bp, while that of B. tuberogranulata was 16,841 bp in length after assembly using NOVOPlasty (Figure 2).Blastn searches revealed that the most similar sequences in NCBI to be those from Boulenophrys and the best-matched sequences showed the species identification was correct.Both mitogenomes exhibited the typical characteristics of animal mitogenomes, comprising 37 genes, including 13 PCGs, 22 tRNAs, two rRNAs, and a D-loop control region (Table 1 1).Nucleotide composition analysis of the 13 PCGs revealed similar negative AT and GC skew patterns, except for the ND6 gene, which exhibited an extraordinary but positive GC skew (Table 3).The D-loop was located between tRNA Pro and tRNA Phe in both species, and it exhibited similar A + T contents, AT skew, and GC skew patterns with PCGs rather than rRNAs and tRNAs (Table 3).

| Positive selection and phylogenetic relationships
The Ka/Ks ratio was calculated to assess the evolutionary rates of each PCG (Figure 4).The highest Ka/Ks value was observed in ATP8 (0.246), while the lowest was in COI (0.037).Other PCGs such as ATP6, ND2, and ND4 exhibited relatively fast evolutionary rates, whereas COI, COIII, and CYTB showed relatively slow rates.
However, all Ka/Ks ratios for the 13 PCGs were less than 1, indicating that they did not exhibit strong positive selection signals and possibly evolved under purifying selection.

| DISCUSS ION
To the best of our knowledge, this study represents the first assembly of the mitogenomes for two horned toads, B. sangzhiensis and B. tuberogranulata.The exploration of similarities and differences in gene orders, genetic structures, base compositions, evolutionary features, and codon usage offers valuable molecular insights into the taxonomic and phylogenetic characteristics of closely related species (Sun et al., 2022).Overall, the mitogenomes of the two species exhibited similarities in terms of size, organization (Figure 2 and Table 1), nucleotide composition of PCGs, rRNAs, tRNAs, control region, and codon usage of PCGs (Figures 3 and 4 and Table 3).The most notable difference between the two mitogenomes was observed in the length of the D-loop (1496 bp vs. 1380 bp), contributing significantly to the total length variations (16,950 bp vs. 16,841 bp).
Despite this distinction, the two newly obtained mitogenomes exhibited patterns similar to those of other Boulenophrys species reported previously (Liu et al., 2016;Wu et al., 2024), such as negative AT skew and GC skew, along with high A + T content.The AT skew and GC skew generally reflect base bias, which can vary among different groups.For instance, the species we previously detected in Salamandridae had a positive AT skew (Wang, Lan et al., 2022), whereas all the Megophryidae species in this study displayed a negative AT skew, consistent with the findings of Zhou et al. (2023).
Furthermore, while the mitogenomes of species in Ranoidea exhibited gene rearrangements frequently (Igawa et al., 2008), no any rearrangements were identified among all the Megophryidae species examined in this study.
By utilizing the mitogenomes of 26 species available from NCBI and the two newly obtained in this study as ingroups, with M. fissipes from the Microhylidae family as the outgroup, we successfully reconstructed a phylogenetic tree within Megophryidae based on the 13 PCGs of 29 species (Figure 5).The phylogenetic analysis revealed that Megophryidae could be divided into two well-supported clades, supporting the classification of the two subfamilies, Megophryinae and Leptobrachiinae.The intergeneric relationships among the two subfamilies were consistent with our previous study (Zhou et al., 2023), with the monophyly of each examined genus also receiving support.
Based on the two most recent phylogenetic studies of Boulenophrys, the genus can be further divided into three well-supported species   (Lyu et al., 2023;Qi et al., 2021) tuberogranulata, B. jinggangensis, B. boettgeri, B. kuatunensis, and B. baishanzuensis grouped in the B. boettgeri group.These phylogenetic relationships within Boulenophrys, as elucidated here using mitogenomes, were largely consistent with previous studies that employed more species but fewer molecular markers (Lyu et al., 2023;Qi et al., 2021).Leptobrachium we studied before, where L. boringii and L. liui were also collected in the same region of Zhangjiajie City but occupied very different phylogenetic lineages that showed a parallel southwestern and southeastern convergence (Zhou et al., 2023).Notably, another species, B. caudoprocta, was also found in Mt.Tianping in the same area (Figure 1).Although the phylogenetic position of B. caudoprocta was not assessed in this study, a previous investigation indicated that the three species found in Mt.Tianping did not form any sister group relationships and they were also morphologically distinct (Lyu et al., 2023).This distribution pattern in Boulenophrys is unique, with three species coexisting in the same mountain but having distinct phylogenetic histories.This pattern highlights the special role of the coexisting area, the Zhangjiajie City in driving the speciation of these species.Similar coexistence of species in nearby areas was also observed in other species of Boulenophrys.were all found in Mt.Jinggang (Jiangxi Province, China).However, these sympatric species were all clustered in the same phylogenetic clade (Lyu et al., 2023).
There are more than 300 species in Megophryidae; however, the available mitogenome data are still very limited.As of now, only 25 complete and nine nearly complete mitogenomes belonging to 26 recognized species of Megophryidae are available from NCBI.The phylogenetic tree presented in this study, based on 28 mitogenomes in Megophryidae, demonstrated that intergeneric and intrageneric relationships can be reliably elucidated with high supporting values.
It suggests that mitogenomes serve as valuable molecular markers for constructing a robust phylogenetic tree of Megophryidae.It is worth noting that the advent of NGS has revolutionized genomics research by enabling the sequencing of entire genomes with everincreasing throughput and decreasing costs (Van Dijk et al., 2014).
With the application of NGS, we anticipate that the complete evolutionary history of Megophryidae will gradually be unveiled when more mitogenomes, such as the two Boulenophrys species in this study, become available in the future.
lections and experiment protocols were approved by the Biomedical Ethics Committee of Jishou University (No: JSDX-2024-0083) adhered to the relevant laws and guidelines of China.Following the "3R principle" (Reduction, Replacement, and Refinement) of animal ethics that required by National Ministry of Science and Technology (No. 398 [2006]), only one sample of each species was utilized.The sample of B. sangzhiensis was collected in May 2021 from its type locality within the National Nature Reserve of Badagongshan in Sangzhi County.The sample of B. tuberogranulata was collected in August 2022, also from its type locality within the Zhangjiajie National Forest Park in Wulingyuan District.The straight-line distance between the two sample localities was only 60 km, and both are situated in Zhangjiajie City, Hunan Province, China.All animal collection and treatment procedures strictly followed the guidance outlined in the Code of Practice for the Housing and Care of Animals.The specimens were euthanized and preserved in 85% alcohol as voucher specimens, and then deposited at the molecular ecology laboratory in Zhangjiajie Campus, Jishou University (B.sangzhiensis, voucher no.JWS20210007; B. tuberogranulata, voucher no.JWS20221040).A small section of liver tissue from each sample was used for molecular analysis.Total DNA extraction was performed using the DNeasy Blood & Tissue Kit (Qiagen), and the DNA library F I G U R E 1 The recorded distribution sites of Boulenophrys species in China (data from AmphibiaChina, 2024).The circles (in red show the species we collected in this study) on the map represent the type localities of Boulenophrys species distributed in China or near the national boundaries.
the VAHTS Universal DNA Library Prep Kit for Illumina V3 (Vazyme).High-throughput sequencing was conducted in paired-end mode on the DNBSEQ-T7 platform (Complete Genomics and MGI Tech), generating approximately 30 Gb of raw reads with a read length of 150 bp for each sample.

3. 3 |
Characteristics of rRNAs, tRNAs, and the control region Both B. sangzhiensis and B. tuberogranulata had two rRNA genes: the 12S rRNA located between tRNA Phe and tRNA Val with a length of 925 bp in both species, and the 16S rRNA located between tR-NA Val and tRNA Leu with lengths of 1594 bp in B. sangzhiensis and 1595 bp in B. tuberogranulata.Both mitogenomes contained 22 tRNAs, with total lengths of 1529 and 1525 bp for B. sangzhiensis, and B. tuberogranulata, respectively.Individual tRNA sizes ranged from 62 bp (tRNA Ser1 ) to 75 bp (tRNA Leu2 ) in B. sangzhiensis and from 63 bp (tRNA Cys ) to 74 bp (tRNA Leu2 ) in B. tuberogranulata.Unlike the pattern of PCGs, both rRNAs and tRNAs had positive AT skew, and the tRNAs additionally exhibited positive GC skew, indicating a different base bias compared to PCGs.The noncoding control region, also known as the D-loop, showed significant size variations, with a length of 1496 bp in B. sangzhiensis and 1380 bp in B. tuberogranulata.

(
Figure 5).The Megophryidae family could be divided into two major clades, corresponding to the subfamilies Leptobrachiinae (Clade I) and Megophryinae (Clade II).Clade I could be further divided into four well-supported major groups: Leptobrachella, Scutiger, Oreolalax, and Leptobrachium.Clade II contained species from the genera Atympanophrys, Brachytarsophrys, and Boulenophrys.Each examined genus appeared as monophyletic, and the intergeneric relationships among the two clades were also well-supported.The genus F I G U R E 2 Gene map of the mitogenome of Boulenophrys sangzhiensis (a) and Boulenophrys tuberogranulata (b).Boulenophrys could be further divided into two subgroups, with the two species sequenced in this study appearing in each of these subgroups.Boulenophrys sangzhiensis grouped with B. omeimontis and then clustered with B. spinata.Comparatively, B. tuberogranulata clustered with B. jinggangensis, B. boettgeri, B. kuatunensis, and B. baishanzuensis.

F I G U R E 3
Relative synonymous codon usage (RSCU) of the mitogenome of Boulenophrys sangzhiensis (a) and Boulenophrys tuberogranulata (b).The Ka/Ks of 13 PCGs among 28 species within Megophryidae.F I G U R E 5 Phylogenetic relationships within Megophryidae derived from BI method based on 13 PCGs.The numbers on the branches represent the bootstrap values and posterior probabilities of ML/BI analyses.The number after species name is the GenBank accession number.Names in red show the phylogenetic positions of Boulenophrys sangzhiensis and Boulenophrys tuberogranulata that we sequenced in this study.
B. minor group, the B. omeimontis group, and the B. boettgeri group Interestingly, despite B. sangzhiensis and B. tuberogranulata being sympatric and coexisting in Zhangjiajie City, Hunan Province, this study revealed that they belong to two distinct phylogenetic clades (Figure 5).Boulenophrys sangzhiensis, along with other species in the B. omeimontis species group, mainly displayed a southwestern distribution pattern.In contrast, B. tuberogranulata and other species in the B. boettgeri species group were primarily distributed in the southeastern area of China.It appears that Zhangjiajie City serves as the boundary area between the two clades of Boulenophrys, especially considering that the type localities of B. sangzhiensis and B. tuberogranulata completely overlapped in Mt.Tianping in Sangzhi County.A similar pattern was observed in mustache toads of For instance, the type localities of B. boettgeri, B. ombrophila and B. kuatunensis were all in Kuatun village (Fujian Province, China), B. yangmingensis and B. xiangnanensis were both in Mt.Yangming (Hunan province, China), and B. cheni, B. jinggangensis, and B. lini Characteristics of the mitogenome of Boulenophrys sangzhiensis (BS) and Boulenophrys tuberogranulata (BT).
Mitochondrial genomes (or called mitogenomes) serve as valuable molecular markers and have found widespread applications in molecular biology and ecology.In vertebrates, especially, mitogenomes are extensively utilized due to their conservative structure, typically TA B L E 1 Base composition (in percentages) of the mitogenomes of Boulenophrys sangzhiensis, Boulenophrys tuberogranulata, and other 26 species in Megophryidae.Nucleotide composition and skewness of different gene regions in the mitogenomes of Boulenophrys sangzhiensis (BS) and Boulenophrys tuberogranulata (BT).
TA B L E 2 aThe sequence obtained in this study.TA B L E 3 . In this study, two species groups were identified and supported, although no sequences of the B. minor group were included.Specifically, B. sangzhiensis, B. omeimontis, and B. spinata formed the B. omeimontis group, while B.